import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import random

# 生成a到b之间的随机数
def GetRandomFactor():
    random.seed()  # 使用系统时间等随机源初始化随机数生成器
    return random.uniform(0.9990, 1.0001)

def GetRandomNumber():
    random.seed()  # 使用系统时间等随机源初始化随机数生成器
    return random.uniform(-0.1, 0.1)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用 SimHei 字体
plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题

# 定义季节函数
def season_factor(t):
    # 将一年分为四个季节，每个季节分别为3个月
    season = (t % 12) // 3  # 每个季节3个月
    return season

# 定义随时间变化的参数
def r(t):
    season = season_factor(t)

    if season == 1:
        increase_rate = 0.2
    elif season == 2:
        increase_rate = 0.24
    elif season == 3:
        increase_rate = 0.25
    else:
        increase_rate = 0.16
    return increase_rate * GetRandomFactor() * 0.7 + 1.2 * np.sin(0.005 * t)

# 定义微分方程组
def evolution_forests_to_farmland(y0, t, params):
    P, I, B, S, AA, BB = y0  # 设置初始种群数量
    K, a, C, a_P, b, I_M, yp_I, rr, yp, Ba, B_M, mu_B, k, S_M, lam, a_AA, m, et_BB, r_AA, K_AA, n, mu_AA, r_BB, K_BB, mu_BB = params

    '''
    K: 植物环境容纳量
    a: 农药对植物影响系数
    C: 化学药剂浓度
    a_P: 植物对昆虫的抗性系数
    b: 昆虫增长率
    I_M: 昆虫种群最大承载量
    yp_I: 昆虫自然死亡率
    rr: 化学物对昆虫影响系数
    yp: 鸟捕虫的效率
    Ba: 鸟类增长率
    B_M: 鸟类最大承载量
    mu_B: 鸟自然死亡率
    k: 土壤恢复速率
    S_M: 土壤最大肥力
    lam: 化学物对土壤影响
    a_AA: 物种AA对植物的食害系数
    m: 物种AA对昆虫的竞争系数
    et_BB: 物种BB对昆虫的捕食系数
    r_AA: 物种AA的增长率
    K_AA: AA的环境容纳量
    n: 昆虫对物种AA的竞争系数
    mu_AA: AA的自然死亡率
    r_BB: BB的自然增长率
    K_BB: BB的环境容纳量
    mu_BB: BB的自然死亡率
    '''
    dPdt = r(t) * P * (1 - P / K) - a * P * C - a_P * P * I - a_AA * P  # 植物增长
    dIdt = b * I * (1 - I / I_M - m * AA / I_M) - rr * I * C - yp * B * I - yp_I * I - et_BB * I + (P / 300) * 15  # 昆虫增长
    dBdt = Ba * B * (1 - B / B_M) + yp * B * I - mu_B * B  # 鸟类增长
    dSdt = k * (S_M - S) - lam * S * C  # 土壤质量
    dAAdt = r_AA * AA * (1 - AA / K_AA - n * I / K_AA) + a_AA * P - mu_AA * AA  # 物种AA增长
    dBBdt = r_BB * BB * (1 - BB / K_BB) + et_BB * I - mu_BB * BB  # 物种BB增长
    return [dPdt, dIdt, dBdt, dSdt, dAAdt, dBBdt]

# 初始种群数量
P0 = 350  # 植物
I0 = 130  # 昆虫
B0 = 32  # 鸟类
S0 = 1.2  # 土地肥力
AA0 = 37  # 物种AA
BB0 = 42  # 物种BB
y0 = [P0, I0, B0, S0, AA0, BB0]

# 时间点
t = np.linspace(0, 150, 1000)

# 参数设计
params = (
    500,  # K: 植物环境容纳量
    2.8,  # a: 农药对植物影响系数
    0.01,  # C: 化学药剂浓度
    0.0005,  # a_P: 植物对昆虫的抗性系数
    0.38,  # b: 昆虫增长率
    200,  # I_M: 昆虫种群最大承载量
    0.05,  # yp_I: 昆虫自然死亡率
    20.8,  # rr: 化学物对昆虫影响系数
    0.0007,  # yp: 鸟捕虫的效率
    0.1,  # Ba: 鸟类增长率
    35,  # B_M: 鸟类最大承载量
    0.08,  # mu_B: 鸟自然死亡率
    0.13,  # k: 土壤恢复速率
    1.5,  # S_M: 土壤最大肥力
    2.8,  # lam: 化学物对土壤影响

    0.015,  # a_AA: 物种AA对植物的食害系数
    0.01,  # m: 物种AA对昆虫的竞争系数
    0.15,  # et_BB: 物种BB对昆虫的捕食系数
    0.2,  # r_AA: 物种AA的增长率
    50,  # K_AA: 物种AA的环境容纳量
    0.02,  # n: 昆虫对物种AA的竞争系数
    0.18,  # mu_AA: 物种AA的自然死亡率
    0.3,  # r_BB: 物种BB的自然增长率
    50,  # K_BB: 物种BB的环境容纳量
    0.3  # mu_BB: 物种BB的自然死亡率
)

# 求解微分方程组
solution = odeint(evolution_forests_to_farmland, y0, t, args=(params,))

# 绘制结果
plt.figure(figsize=(10, 8))
plt.plot(t, solution[:, 0] / P0, label='植物当前生物量/初始生物量 (P)')
plt.plot(t, (solution[:, 1] / I0) * GetRandomFactor() + GetRandomNumber() * 0.55 * GetRandomFactor() * 1.2*np.sin(t), label='昆虫当前生物量/初始生物量 (I)')
plt.plot(t, solution[:, 2] / B0 + GetRandomNumber() * 0.05 * GetRandomFactor() * np.sin(t), label='原生鸟类当前生物量/初始生物量(B)')
plt.plot(t, solution[:, 3] / S0, label='土地肥力 /初始土地肥力(S)')
plt.plot(t, solution[:, 4] / AA0+ GetRandomNumber() * 0.05 * GetRandomFactor() * np.cos(t), label='食草性动物A当前生物量 /初始生物量(A)')
plt.plot(t, solution[:, 5] / BB0+ GetRandomNumber() * 0.05 * GetRandomFactor() * np.sin(t), label='鸟类B当前生物量 /初始生物量(B)')

plt.legend()

plt.xlabel('时间（月）')
plt.ylabel('初始生物量（土地肥力）/当前生物量（土地肥力）')
plt.title('农田生态系统成熟后种群规模/土地肥力随时间变化（使用化学物）')
plt.show()

